# include "scene.h"
# include <stdio.h>
# include <math.h>
# include <stdlib.h>

#ifdef _OMP_ORLANSKY_OUTLET_
# include <omp.h>
#endif

struct back_flow{ unsigned short i,j; char flag; };

extern double orlansky_outlet(struct grid *g, unsigned int t_cnt, const struct arr_struct *memOut)
   {
    short i, j, k;
    unsigned short tn, nr, nth, nTRD, tid, i_min, j_min;
    double cdtdz, dt, *th, *r, *z, q0, q1, m_inlet, m_accum, ***u, ***v, ***w;
#   ifdef SOLVE_ENERGY
    double ***T;
#   endif
    struct back_flow *bf;
    char bf_flag;

    r = g->r[2]; nr = g->nr[2]; nth = g->nth; th = g->th; z = g->z[2]; dt = g->dt;
    u = g->u[2]; v = g->v[2]; w = g->w[2];
#   ifdef SOLVE_ENERGY
    T = g->T[2];
#   endif

    nTRD = 1;
#   ifdef _OMP_ORLANSKY_OUTLET_
    nTRD = g->nTRD;
#   endif

 /* Allocating memory for the structure that reflects if back-flow has taken place or not. */
    bf = (struct back_flow *)malloc(sizeof(struct back_flow)*nTRD);

 /* Mass coming into the domain through inlet between time t~t+dt. */
    m_inlet = 0.25*g->density0*r_NOZZLE*r_NOZZLE*W_MAX\
       *(dt-tFRAC*(exp(log(1-wFRAC)*(g->t+dt)/tFRAC) - exp(log(1-wFRAC)*g->t/tFRAC))/log(1-wFRAC))*(th[NthG+g->nth]-th[NthG]);

 /* Rearraning the memory to backup the physical quantities at current time-instant. */
    for(i = t_cnt%tn_BACKUP, tn = k = 0; tn < tn_BACKUP; tn++)
         {
          if(tn == i) k = tn_BACKUP;
          g->outBC.u[tn] = memOut->q[  tn_BACKUP+tn-i-k];
          g->outBC.v[tn] = memOut->q[2*tn_BACKUP+tn-i-k];
          g->outBC.w[tn] = memOut->q[3*tn_BACKUP+tn-i-k];
#         ifdef SOLVE_ENERGY
          g->outBC.T[tn] = memOut->q[4*tn_BACKUP+tn-i-k];
#         endif
         }

  q0 = (double)0;
# ifdef _OMP_ORLANSKY_OUTLET_
  #pragma omp parallel private(i,j,k,cdtdz,q1,tid)
# endif
   {
    tid = 0;
#   ifdef _OMP_ORLANSKY_OUTLET_
    tid = omp_get_thread_num();
#   endif
    bf[tid].flag = 'n';

 /* Backing up the quantities at current time-instant. */
#   ifdef _OMP_ORLANSKY_OUTLET_
    #pragma omp for schedule(guided)
#   endif
    for(i = 0; i < nr; i++) for(k = 0; k < zCells_BACKUP; k++) for(j = 0; j < nth; j++)
       {
        g->outBC.u[0][i][j][k] = u[NrG+i][NthG+j][NzG+k-1];
        g->outBC.v[0][i][j][k] = v[NrG+i][NthG+j][NzG+k-1];
        g->outBC.w[0][i][j][k] = w[NrG+i][NthG+j][NzG+k-1];
#       ifdef SOLVE_ENERGY
        g->outBC.T[0][i][j][k] = T[NrG+i][NthG+j][NzG+k-1];
#       endif
       }

 /* Calculating quantities at the first adjacent ghost cells near the outlet at next time step using the Orlansky BC. */
#   ifdef _OMP_ORLANSKY_OUTLET_
    #pragma omp for schedule(guided)
#   endif
    for(i = 0; i < nr; i++) for(j = 0; j < nth; j++)
       {
     /* Calculating "tnext_ug[i][j]" at the outlet. */
        cdtdz = g->outBC.u[2][i][j][1] - g->outBC.u[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.u[0][i][j][1] + g->outBC.u[2][i][j][1] - 2*g->outBC.u[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_ug[i][j] = (1-cdtdz)*g->outBC.u[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.u[0][i][j][1]/(1+cdtdz);

     /* Calculating "tnext_vg[i][j]" at the outlet. */
        cdtdz = g->outBC.v[2][i][j][1] - g->outBC.v[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.v[0][i][j][1] + g->outBC.v[2][i][j][1] - 2*g->outBC.v[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_vg[i][j] = (1-cdtdz)*g->outBC.v[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.v[0][i][j][1]/(1+cdtdz);

     /* Calculating "tnext_wg[i][j]" at the outlet. */
        cdtdz = g->outBC.w[2][i][j][1] - g->outBC.w[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.w[0][i][j][1] + g->outBC.w[2][i][j][1] - 2*g->outBC.w[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_wg[i][j] = (1-cdtdz)*g->outBC.w[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.w[0][i][j][1]/(1+cdtdz);
        if(bf[tid].flag != 'y') if(g->outBC.tnext_wg[i][j] > 0)
            {
             bf[tid].i = NrG+i;
             bf[tid].j = NthG+j;
             bf[tid].flag = 'y';
            }

     /* Calculating "tnext_Tg[i][j]" at the outlet. */
#       ifdef SOLVE_ENERGY
        cdtdz = g->outBC.T[2][i][j][1] - g->outBC.T[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.T[0][i][j][1] + g->outBC.T[2][i][j][1] - 2*g->outBC.T[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_Tg[i][j] = (1-cdtdz)*g->outBC.T[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.T[0][i][j][1]/(1+cdtdz);

// fprintf(g->lptr, "%d, %d, %E\n", i, j, g->outBC.tnext_Tg[i][j]);
#       endif
       }
// if(t_cnt == 1) exit(-1);

 /* Mass going out of the domain between time t~t+dt. */
#   ifdef _OMP_ORLANSKY_OUTLET_
    #pragma omp for reduction(+:q0)
#   endif
    for(i = 0; i < nr; i++)
        {
         for(q1 = 0, j = 0; j < nth; j++)
             q1 += (th[NthG+j+1]-th[NthG+j])*(g->outBC.w[0][i][j][0]+g->outBC.tnext_wg[i][j]);

         q0 += (r[NrG+i+1]+r[NrG+i])*(r[NrG+i+1]-r[NrG+i])*q1;
        }

#   ifdef _OMP_ORLANSKY_OUTLET_
    #pragma omp single
#   endif
      {
    /* m_accum = m_inlet - m_outlet. */
       m_accum = m_inlet - (-0.25*g->density0*q0*dt);

    /* Determining the correction factor for "w_outlet" so that "m_accum = 0" between time t~t+dt. */
       q0 = 0;
#      ifdef N_w_CORRECT
       if(t_cnt < N_w_CORRECT)
#      endif
       q0 = 4*m_accum/(g->density0*dt*(r[NrG+nr]*r[NrG+nr]-r[NrG]*r[NrG])*(th[NthG+nth]-th[NthG]));
      }

 /* Correcting the w_outlet.  */
#   ifdef _OMP_ORLANSKY_OUTLET_
    #pragma omp for schedule(guided)
#   endif
    for(i = 0; i < nr; i++) for(j = 0; j < nth; j++) g->outBC.tnext_wg[i][j] -= q0;
   }

 /* Checking if the back-flow has occurred or not. */
    for(bf_flag = 'n', tid = 0; tid < nTRD; tid++) if(bf[tid].flag == 'y')
       {
        if(bf_flag == 'y')
            {
             if(bf[tid].i < i_min)
                 {
                  i_min = bf[tid].i;
                  j_min = bf[tid].j;
                 }
             else if(bf[tid].i == i_min) if(bf[tid].j < j_min) j_min = bf[tid].j;
            }
        else{
             i_min = bf[tid].i;
             j_min = bf[tid].j;
             bf_flag = 'y';
            }
       }

    if(bf_flag == 'y')
       {
        fprintf(g->lptr, "ERROR:: In 'orlansky_outlet()': Backflow observed at outlet.\n\tLocation: (i,j) = (%d, %d)\n\tLocation: (r,th) = (%E mm, %E Degrees)\n\tTime-instant: t = %E sec;\n\tt_cnt = %d\nExiting.\n", i_min, j_min, 0.5*(r[i_min]+r[i_min+1])*1000, 45*(th[j_min]+th[j_min+1])/acos(0), g->t+dt, t_cnt+1);
        exit(-1);
       }

    free(bf);

 /* Return % of dmdt-accumulation at next time-stage. */
    return m_accum*100/m_inlet;
   }
